// rewrite from https://oi-wiki.org/math/number-theory/pollard-rho/

use rand::Rng;

pub fn find_max_prime_factor(number: u128) -> u128 {
    let mut max_factor = 0u128;
    fac(number, &mut max_factor);
    return max_factor;
}

fn gcd(mut a: u128, mut b: u128) -> u128 {
    while b != 0 {
        let t = a % b;
        a = b;
        b = t;
    }
    a
}

fn bmul(a: u128, b: u128, m: u128) -> u128 {
    // an overflow version
    // (a * b) % m
    let mut ans = 0u128;
    let mut a = a % m;
    let mut b = b % m;
    while b > 0 {
        if b % 2 == 1 {
            ans = (ans + a) % m;
        }
        a = (a + a) % m;
        b /= 2;
    }
    ans
}

fn qpow(mut x: u128, mut p: u128, m: u128) -> u128 {
    let mut ans = 1u128;
    x %= m;
    while p > 0 {
        if p % 2 == 1 {
            ans = bmul(ans, x, m);
        }
        x = bmul(x, x, m);
        p /= 2;
    }
    ans
}

fn miller_rabin(p: u128) -> bool {
    if p < 2 {
        return false;
    }
    if p == 2 || p == 3 {
        return true;
    }
    if p % 2 == 0 {
        return false;
    }
    let mut d = p - 1;
    let mut r = 0;
    while d % 2 == 0 {
        d /= 2;
        r += 1;
    }
    let mut rng = rand::thread_rng();
    for _ in 0..10 {
        let a = rng.gen_range(2..p - 1);
        let mut x = qpow(a, d, p);
        if x == 1 || x == p - 1 {
            continue;
        }
        let mut composite = true;
        for _ in 0..r - 1 {
            x = bmul(x, x, p);
            if x == p - 1 {
                composite = false;
                break;
            }
        }
        if composite {
            return false;
        }
    }
    true
}

fn pollard_rho(x: u128) -> u128 {
    if x % 2 == 0 {
        return 2;
    }
    let mut rng = rand::thread_rng();
    let c = rng.gen_range(1..x);
    let mut s = 0u128;
    let mut t = 0u128;
    let mut val = 1u128;
    let mut step = 0;
    let mut goal = 1;
    loop {
        for _ in 0..goal {
            t = (bmul(t, t, x) + c) % x;
            val = bmul(val, t.abs_diff(s), x);
            step += 1;
            if step % 127 == 0 {
                let d = gcd(val, x);
                if d > 1 && d < x {
                    return d;
                } else if d == x {
                    return pollard_rho(x);
                }
            }
        }
        let d = gcd(val, x);
        if d > 1 && d < x {
            return d;
        } else if d == x {
            return pollard_rho(x);
        }
        s = t;
        goal *= 2;
        val = 1;
    }
}

fn fac(x: u128, max_factor: &mut u128) {
    if x <= *max_factor || x < 2 {
        return;
    }
    if miller_rabin(x) {
        *max_factor = std::cmp::max(*max_factor, x);
        return;
    }
    let mut p = x;
    while p >= x {
        p = pollard_rho(x);
    }
    let mut x = x;
    while x % p == 0 {
        x /= p;
    }
    fac(x, max_factor);
    fac(p, max_factor);
}
